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We introduce a method to obtain the specific heat of quantum impurity models via a direct 
calculation of the impurity internal energy requiring only the evaluation of local quantities within a 
single numerical renormalization group (NRG) calculation for the total system. For the Anderson 
impurity model, we show that the impurity internal energy can be expressed as a sum of purely 
local static correlation functions and a term that involves also the impurity Green function. The 
temperature dependence of the latter can be neglected in many cases, thereby allowing the impurity 
specific heat, Ci mp , to be calculated accurately from local static correlation functions; specifically 

via Ci m p = dE i j ?p"' + \ 9 g'y b , where Skmic and -Bhyb are the energies of the (embedded) impurity and 
the hybridization energy, respectively. The term involving the Green function can also be evaluated 
in cases where its temperature dependence is non- negligible, adding an extra term to Ci mp . For 
the non-degenerate Anderson impurity model, we show by comparison with exact Bethe ansatz 
calculations that the results recover accurately both the Kondo induced peak in the specific heat 
at low temperatures as well as the high temperature peak due to the resonant level. The approach 
applies to multiorbital and multichannel Anderson impurity models with arbitrary local Coulomb 
interactions. An application to the Ohmic two state system and the anisotropic Kondo model is 
also given, with comparisons to Bethe ansatz calculations. The approach could also be of interest 
within other impurity solvers, for example, within quantum Monte Carlo techniques. 

PACS numbers: 75.20.Hr, 71.27.+a, 72.15.Qm 



I. INTRODUCTION 

Quantum impurity models play an important role in 
condensed matter physics, for example, as models of tran- 
sition metal and rare-earth impurities in metals or two- 
level systems and qubits interacting with an envi- 
ronment or in describing the Kondo effect in nanoscale 
devices such as molecular transistors, 8-11 semiconductor 
quantum dots, carbon nanotubes, 15 and magnetic 
ions such as Co 16 ' 11 or Ce ls adsorbed on surfaces. In ad- 
dition, they appear as the effective models within dynam- 
ical mean field theory (DMFT) treatments of strongly 
correlated electron systems, such as heavy fermions and 
transition metal oxides. 19-22 Hence, new approaches to 
calculate their dynamic, thermodynamic and transport 
properties are potentially of wide interest. 

The numerical renormalization group (NRG) 
method, in particular, has proven very suc- 

cessful for the study of quantum impurity models. The 
method, described briefly in the next section, gives both 
the thermodynamic, 23_25 > 2 ' dynamic, 28-35 and transport 
properties '' 1 of quantum impurities. Thermodynamic 
properties, such as the specific heat, are of particular 
interest for bulk systems, such as dilute concentrations 
of transition metal or rare-earth ions in non-magnetic 
metals. A measurement of the temperature dependence 
of the specific heat or susceptibility of such systems 
provides important information about their physical 
behavior, for example, whether such systems exhibit 
Fermi liquid or non-Fermi liquid behavior at low tem- 
perature and thus information about the nature of their 
low energy excitations. ' 



The usual approach to calculating the specific heat of 
quantum impurity models within the NRG method con- 
sists of a two-stage procedure, 24-2 ' in which the Hamil- 
tonians of the total system H is first diagonalized, fol- 
lowed by a similar diagonalization for the host Hamilto- 
nian H . Here, H = H smp + H int + H is the Hamil- 
tonian of a quantum impurity (described by H- lmp ), in- 
teracting with a host (described by H ) via the interac- 
tion term H lnt - From the eigenvalues of H and Hq, the 
grand canonical partition functions Z = Tr e~@ H and 
Zq = Tr e~P H ° and the corresponding thermodynamic 
potentials Q(T) = -k B T\nZ and Q (T) = -k B T\nZ a 
are constructed, where /3 — l/k B T is the inverse tem- 
perature. The impurity contribution to the specific heat, 
Ci mp {T), is then obtained by subtraction via Ci mp (T) = 
C(T) - C (T), where C(T) and C (T) are the specific 
heats of the total system and of the host system, respec- 
tively, 

C ( T ) = - T& ^j? = ^((H < H)f), (1) 

C o(T) = -T &2 f^ = k B f((H Q -(H )Y) (2) 

C lmp (T) = C(T) - C (T). (3) 

In this paper we present a new approach to the calcula- 
tion of the impurity internal energy and specific heat of 
quantum impurity models within the numerical renor- 
malization group (NRG) method. 23-26 It relies on ex- 
pressing the impurity internal energy in terms of local 
quantities, and as such is not restricted to the NRG but 
may be implemented within any impurity solver that cal- 
culates such quantities. The main result of this paper is 
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the (approximate) expression for the impurity specific 
heat of the Anderson model (see Sec. Ill) 
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vantages of this approach are that, (i), Eq. (4) involves 
only a first temperature derivative and is expected to be 
more accurate for numerical evaluations than Eqs. (1)- 
(3) which involve a second temperature derivative of the 
thermodynamic potential, or, the calculation of the total 
energy fluctuation, (ii), the host contribution to the in- 
ternal energy (Ho) has been analytically subtracted out 
(see Sec. Ill), so only the diagonalization of H is required, 
(iii), only local static correlation functions appearing in 
(Hi mp ) and (H- m t) are required, and, (iv), as we shall 
show, the new approach is less sensitive to discretiza- 
tion effects of the host than the usual approach which 
evaluates expectation values of extensive quantities. We 
illustrate the method by applying it to the Anderson im- 
purity model and we compare the results for specific heats 
with those from the conventional NRG approach 2 '' 36 ' 39 
and with exact results from thermodynamic Bethe ansatz 
calculations . 41 1-4 2 

Early approaches to the specific heat of dilute Kondo 
systems used an equation of motion decoupling scheme 
for the Kondo model and expressed the impurity in- 
ternal energy in terms of the local T-matrix. The re- 
sults obtained for the specific heat within this approx- 
imation were inadequate, violating, for example, Fermi 
liquid properties at low temperatures. 41 A formally ex- 
act expression for the internal energy of the Anderson 
model, in terms of the local self-energy and the local 
Green function, was obtained by Kjollerstrom et. al., 
in Ref. 45. They evaluated the specific heat in the low 
density limit (corresponding to a small occupation of the 
local level) obtaining correct results obeying Fermi liquid 
theory in this limit. 

The most reliable approaches to specific heats of quan- 
tum impurity models are the Bethe ansatz method for 
integrable models 40-42 ' 46-49 and the NRG method. An 
important aspect of the latter, allowing it to access ther- 
modynamic properties on all temperature scales down 
to T — 0, is the use of a logarithmic grid to represent 
the quasi-continuous spectrum u e [— D, +D] of the host 
system, H . Thus ui — > ui n = ±DA~ n ,n = 0,1,..., 
where the parameter A > 1 achieves a separation of the 
many energy scales in Ho and thus in H (see Sec. II). 
A large A 3> 1 allows calculations to reach low temper- 
atures in fewer steps within the iterative diagonalization 
procedure of the NRG, and, in addition, a large A ^> 1 
reduces the size of the truncation errors at each step in 
this procedures ! However, for A ^ 1, specific heats (and 
also susceptibilities) , calculated by using a standard loga- 
rithmic grid, exhibit discretization oscillations, especially 
at low temperatures. '" On the other hand, calculations 
at smaller A < 3, with less severe discretization oscil- 
lations, are more prone to truncation errors. In order 



to be able to carry out accurate calculations at all tem- 
peratures, using A ^> 1, an averaging over several dis- 
cretizations of the host degrees of freedom has been in- 
troduced which essentially allows exact calculations to be 
carried out. 50 ' 51 With this refinement, the NRG approach 
has been used extensively in calculations of specific heats 
of quantum impurity models, with applications to the 
two-impurity Kondo model 52 ' 53 and the two-channel An- 
derson models. 54 

The paper is organized as follows. In Sec. II, the An- 
derson impurity model is described, and the NRG is out- 
lined together with a brief description of how thermo- 
dynamic properties are conventionally calculated within 
NRG (at A 1). In Sec. Ill, we describe our new ap- 
proach to specific heats of quantum impurity models, us- 
ing the Anderson impurity model as an example (with 
some further details given in Appendix A). The avail- 
ability of exact Bethe ansatz results for this model, 
allows a detailed evaluation of the accuracy of our new 
approach to specific heats. Results at zero and finite 
magnetic fields are presented in Sec. IV for the symmet- 
ric Anderson model. These are compared to both exact 
Bethe ansatz results and results obtained in the conven- 
tional NRG approach. Sec. V contains results for the 
asymmetric model with comparisons to corresponding 
Bethe ansatz calculations. The thermodynamic Bethe 
ansatz (TBA) equations for the Anderson impurity model 
and the details of their numerical solution can be found in 
Appendix B. In Sec. VI we present the generalization to 
multichannel and multiorbital Anderson impurity models 
and to dissipative two state systems. For the Ohmic case, 
results for specific heats are compared to corresponding 
Bethe ansatz results for the equivalent anisotropic Kondo 
model (AKM). Section VII summarizes the main results 
of this paper and discusses possible future applications. 



II. MODEL, METHOD AND CONVENTIONAL 
APPROACH TO THERMODYNAMICS 

We consider the Anderson impurity model, described 
by the Hamiltonian 



H = H, 



Hq + Hij 



The first term, H imp = e dd\d a + Un d -\n d ^, de- 
scribes the impurity with local level energy Sd and on- 
site Coulomb repulsion U, the second term, Ho = 
Sfca e k c ka c kcr> 1S ^ ne kinetic energy of non-interacting 
conduction electrons with dispersion and, the last 

term, H int = J2ka V k( c k<r d <r + d l c k<y), is the hybridiza- 
tion between the local level and the conduction electron 
states, with Vk being the hybridization matrix element. 
We shall also consider the effect of a magnetic field of 
strength B by adding a term Hb = —9HbB S z to H 
where S z is the z-component of the total spin (i.e., im- 
purity plus conduction electron spin), g is the electron 
(/-factor, and /ib is the Bohr magneton. We choose units 
such that g = fiR = 1. 
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The NRG procedure consists of the following steps. 
First, the conduction electron energies — D < £k < D, 
where D is the half-bandwidth, are logarithmically dis- 
cretized about the Fermi level Ef = 0, that is, tu — > e n — 
±DA~ n , n = 0, 1, . . . where A > 1 is a momentum rescal- 
ing factor. We shall also consider generalized discretiza- 
tions denned by a parameter z, such that e — and 
e n = ±DA~ n ~^~ z ' ,n = I,... , with 2 = 1 recovering the 
usual discretization. For A > 1, discretization induced 
oscillations of period In A can be eliminated by averag- 
ing results for several z in (0, 1]. 50 > 51 Second, the opera- 
tors c na , n — 0,1, ... , are rotated to a new set f na , n = 
0, 1, . . . , with Vfoa- = J2^=o v k n c na , such that the dis- 



E„=0<7 ± E n{z)a 



cretized conduction band Hq 

with, for example, E n (z) = ~(1 + A~ L )Dk~ n for z = 1, 

takes the tri-diagonal form H -t J2™=o<j ^nW/L/w + 

E~=o<t *n(f)(/L/n+ia + fl+iafno) in the new basis. Fi- 
nally, within this new basis, the sequence of truncated 
Hamiltonians H m , m — 0, 1, ... , where H m = H lmv + 

-ffhyb + ^2 n =0cr ^ n ( Z )fnafna + Eh=0j ( Z ) (fnafn+lcr + 
fl+lafnv), With H hyh = Vj2a(fL d ° + <4/o<t), IS 

iteratively diagonalized by using the recursion re- 
lation ff m+ l = if m + Yscr em+l( z )/i+l (J /m+l<T + 

E^mWf/L/m+b + fL+tafma)- This procedure 24 -' 2 ' 
yields the eigenstates |p) m and eigenvalues -E™ on a de- 
creasing set of energy scales uj m (z) ~ t m (z), m = 0, 1, 

Since the number of states increases as 4 m+2 , only the 
lowest states are retained for m > m , where typically 
mo > 4 — 5. This is implemented either by, (i), specify- 
ing an approximately constant number of states Nk eep to 
retain at each m > mo, and mo will be fixed by the pre- 
cise value of N-^ccp, or, (ii), by specifying that only those 
states with rescaled energies (E™ — E^ s )/t m (z) < e c (A) 
be retained for m > mo, for some predefined mo, where 
Eqs is the (absolute) groundstate energy at iteration 
m and e c (A) is A-dependent cut-off energy. Combin- 
ing the information from all iterations then allows the 
calculation of thermodynamics on all temperature scales 
of interest.' 39 ' 5 " For most of the results in this paper, 
we used the truncation scheme (ii) with m = 4 — 5 
and e c (A) = 20v / A, similar to the choice in Ref. 39. 
Some calculations using the truncation scheme (i) with 
A?kccp = $60 were also carried out in Sec. VI B. Both 
schemes were found to work well by comparison with ex- 
act Bethe ansatz calculations. Whereas in scheme (i), a 
fixed number, JVk ccp , of levels is retained for all iterations 
m > mo, in scheme (ii), the number of retained states, 
initially large for m < m (typically several thousand), 
starts to decrease with increasing m, eventually saturat- 
ing to a few hundred states at m > mo (e.g., for A = 4). 
While in both schemes only the retained states of itera- 
tion m are used to set up the Hamiltonian H m+ i for the 
next iteration, all states of iteration m are available, and 
are used, in practice, to calculate the thermodynamics. 

The specific heat is calculated within the approach of 
Campo and Oliveira in Ref. - r >l , which we shall refer to 
as the "conventional" approach: For any temperature T, 



| 1.5 
^ 1 

m 

"| 0.5 

c/T „ 



„ 0.2 

m 

1 

°~ -0.2 



-0.4 



""1 1 1 




1 1 iMJ-4-niiii| ill."" 

/ (a) 




— z-averaaed 
-- Z= 1 


- 

'N / 




— ' 


m.l M . 


*~ J I — I 
1 1 1 1 


_ 

1 ' ' ' ""1 1 ' ' ' ""1 1 ' ' ' ""1 1 ' ' ' ""1 1 ' 1 1 llll 


1 11 
1 1 1 
1 1 1 
1 I 1 
1 1 


,1 M 
II " l> 11 

i i , , i 1 ' L 

4-^i ' 1 1 


f\ (b) 


I 1 
- 1 I 
1 I 


' 1 " 
1 i i 1 • 

1 i i 1 


j v \/ 


- 1 1 
1 1 

- M 

M..1 , , 


i i , >i 
i ' i, J 

\l V 

1 1 ' 





10" 



10 1 



io-' 

k B T/D 



10' 



10" 



FIG. 1. (Color online) Temperature dependence of, (a), 
the impurity entropy, S , i mp (T), and, (b), the impurity spe- 
cific heat, Ci m p(T), for the symmetric Anderson model with 
U/A Q = 12 and A = 0.001D. The calculations are for A = 4 
with an energy cut-off e c (A = 4) = 40, without 2-averaging 
[n z = 1, z = 1 (dashed lines)], and with ^-averaging [n z = 2, 
z — 1/4, 3/4 (solid lines)]. For A = 4 two z values suffice to 
eliminate the discretization oscillations. 



we choose the smallest m such that k^T > t m (z) and 
we use the eigenvalues of H m to evaluate the partition 
function Z m {T) = J2 P e~ E ™ /kl>T . The expectation value 
(H) is then calculated, followed by ((H — (H)) 2 ) and the 
specific heat C(T) (in addition, the thermodynamic po- 
tential O(T) = — ksT In Z m (T) may also be calculated). 
Calculations are carried out for several values of the z pa- 
rameter and then averaged. In the calculations reported 
below, we choose z = (2i — \)/2n z ,i = l,...,n z with 
n z = 2,4 or 8. This procedure is repeated for the con- 
duction band Hamiltonian Hq to obtain the host contri- 
bution to the specific heat, Cq(T). Finally, the impurity 
specific heat is obtained via Ci mp (T) = C{T) — Cq(T). 
The above prescription works well for A > 4, since the 
use of large A reduces the size of truncation errors dur- 
ing the iterative diagonalization of H and Hq. m Further- 
more, the use of large A, implies that the highest states 
of H m have energies ^> t m (z) ~ T so that Z m (T) is a 
good approximation to the partition function of the infi- 
nite system at temperature T. In addition to the specific 
heat, we also calculate the impurity contribution to the 
entropy, S imp (T) = S(T)-S (T), where S(T) and S (T) 
are the entropies for H and Hq, respectively, and 

S(T) = -—=k B \nZ(T) + (H)/T, (5) 
an 

So (T) = - ^ = fee In Z Q (T) + (H ) /T. (6) 

Unless otherwise specified, the NRG calculations pre- 
sented in this paper will be for a band of half-width 
D = 1 and a constant particle-hole symmetric density 
of states Np = 1/2D. The hybridization strength, A , 
defined as the half-width of the resonant level is given 
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FIG. 2. (Color online) Temperature dependence of, (a), 
the impurity entropy, 5 , i mp (7 1 ), and, (b), the impurity spe- 
cific heat, Ci mp (T), for the symmetric Anderson model with 
U/A Q = 12 and A = 0.00113. Symbols: NRG calcula- 
tions using the conventional approach. Solid lines: Bethe 
ansatz calculations. The NRG calculations are z-averaged 
with n z = 2 and other parameters as in Fig. 1. 



by Ao = irNpV 2 . Calculations for the positive and 
negative-?/ Anderson models include a U(l) symmetry 
for total electron number conservation and SU(2) sym- 
metry for total spin conservation. We use the discretiza- 
tion scheme of Campo and Oliveira in Ref. 51. 



Figure 1 shows the temperature dependence of the spe- 
cific heat and entropy, calculated with the above proce- 
dure, for the symmetric Anderson model with U/ Aq = 12 
and Ao = 0.001D. The calculations are for A = 4 us- 
ing an energy cut-off e c (A = 4) = 40, both without 
z-averaging (n z — 1) and with z-averaging (n z = 2). 
Note the aforementioned oscillations in the case of no 
z-averaging (n z — 1). For A = 4, two z values suffice 
to eliminate the discretization oscillations (whereas for 
A = 10, four values are required). In order to quan- 
tify the accuracy of the NRG calculations, we also solved 
numerically the thermodynamic Bethe ansatz equations 
for the Anderson model and calculated the entropy and 
specific heat (see Appendix B for details) . A comparison 
of the z-averaged NRG calculations with the exact Bethe 
ansatz results, shown in Fig. 2, indicates very good agree- 
ment. Nevertheless, in the next section we show that the 
specific heat can be calculated directly from the impu- 
rity contribution to the internal energy in terms of lo- 
cal static correlation functions and that discretization ef- 
fects within this approach are less pronounced than those 
above. 



III. IMPURITY INTERNAL ENERGY AND 
SPECIFIC HEATS 

The impurity internal energy is defined by E- lmp 
£totai - E where E total = (H) and E = (H ) 



ka efe\<4r C *W0 



where the subscript denotes a ther- 



modynamic average for non-interacting conduction elec- 
trons (i.e., impurity is absent). We have 



E ° = E / def(e)eN(e), 



(7) 



where /(e) is the Fermi function and N(e) = ^2 k S(e — 
e k ) is the non-interacting conduction electron density of 
states per spin. -Etotai has four contributions: 
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total 



E, 
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E, 



cond 



E 



hyb, 



(8) 



where E occ = J2a £ d( n da), E docc = U(n dt n di ), E cond = 
J2ka e k(cl a c ka ) and E hyh = V J2ka( c L d <? + d\c ka ). The 
first two contributions are evaluated as thermodynamic 
averages within the NRG calculation, requiring the cal- 
culation of matrix elements of ^ CT rida and the double 
occupancy operator D occ = n^ndi- The contribution 
i?i iy b may also be evaluated as a thermodynamic average 
E\yh = VJ2a(dtfocr + H.c). For the discussion below it 
is useful to note that the contribution i?hyb can also be 
expressed in terms of the local retarded d-electron Green 
function Gdo-(w) = {{d a ; dj.))^^ and the hybridization 
function A(u>) = J2k V 2 1 + i5 — e k ) as 



hyb 



-~E/ &jf(u>)Im[G dlT (uj)A 



M] ■ (9) 



Next, 



consider the contribution E r 



This is not simply Eq since the 



impurity affects the conduction electrons once V is 
finite. It can be evaluated from the equation of motion 
of the retarded conduction electron Greens function 

Gfar(w) = {{ckcx\Vk<,))u+iS' 

Gka = G ka + G\ a T a G\ a . (10) 

Here, T a {u) — V 2 Gda(u) is the local T-matrix and 
G ka (ui) — 1/(uj + i5 — e k ) is the non-interacting con- 
duction electron Greens function. Using 



r C ka) 



dujf(uj)lm ((c fctT ;4 CT )) 



we find for E, 



cond 



E, 



cond — E + E lnt 



where 



-Bint = ~~yi J J delm 



eV 2 N(e) 



(u + iS - e) 2 



G dcr (w) 
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where I(co) is given by 
1 



I(lo) 



de- 



cAi(e) 



'(w + i<5 - e) 2 



<9w 



( W A( W )) 



with Ai(e) = Im [A(e + iS)} = -nV 2 N(e), and we evalu- 
ated I(cj) analytically by noting that A(cj + i6) has the 
same properties as a retarded Green function (see Ap- 
pendix A for details). We therefore find, 
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From this and Eq. (9) we see that E- 
Hence, the impurity contribution to the internal energy, 
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-E-docc is adiabatically 



-E 



hyb 



where E\, 



= (H: 



imp/ 



— E c 



E 



connected to the energy of the impurity decoupled from 
the band (i.e., its energy at V — > 0). All contributions 
to E- Imp , except for the last one, can be evaluated as 
thermodynamic averages of local static correlation func- 
tions: The contribution E-Q from the band which in- 
volves a finite frequency Greens function has been related 
to -Ehyb, which can be evaluated as local static correla- 
tion function V J2 a (^tfoa+H.c.) . The contribution E^ , 
also involves a finite frequency Greens function, but we 
could not express this as a local static correlation func- 
tion. Its temperature dependence, however, is negligi- 
ble since the main temperature dependence arises from 
the Fermi window \uj\ < T, but this region is cut out 

(2) 

in i? int due to the factor of w. In addition, for many 
cases of interest d (A(u>)) /doj is small and vanishes in 
the wide band limit: D — > oo and Ao = irN(0)V 2 fixed. 
For example, for a constant density of states it equals 
^(1 - {uj/Dfy 1 ~ A /D for uj < D. Thus, to a 
very good approximation, which we shall quantify in the 
rest of the paper with detailed numerical calculations and 
comparisons to exact Bethe ansatz results, we can ap- 
proximate the impurity contribution to the specific heat 
and entropy via E imp = E ionic + ^E hyb as 
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FIG. 3. (Color online) Temperature dependence of, (a), 
the impurity entropy, Simp (TO, and, (b), the impurity spe- 
cific heat, Ci m p(T), for the symmetric Anderson model with 
U/A = 12 and A = 0.001D calculated within NRG us- 
ing the new approach for A = 4 with an energy cut-off 
e c (A = 4) = 40. Solid lines: n z — 2 (^-averaging). Dashed 
lines: n z = 1 (no z-averaging). For A = 4 two z values thus 
suffice to eliminate the discretization oscillations at n z = 1. 



The omitted term, dEr>/dT, in (16) as argued above, 
has a negligible temperature dependence (although its 
magnitude is not necessarily always small compared to 
the terms retained). Notice that E lmp is made up of 
a term due to the partial occupation of the local reso- 
nant level (E occ ), a term due to the Coulomb repulsion 
of electrons in this level (-Edocc), and, a term due to the 
energy gained by hybridization of the local level with the 
conduction electrons (Ehyb/2), that is, it involves only 
local static correlation functions. Such quantities can 
be calculated very accurately and efficiently within the 
NRG method, within a single calculation for the total 
system only, a significant advantage of this approach. In 
some situations, the hybridization function A(w) may be 
strongly asymmetric and have a strong energy depen- 

(2) 

dence close to u = 0. In such cases, the term E> n £ can be 
calculated via the local spectral function and included 
in -Eimp, which is possible within the NRG, at some- 
what higher numerical cost. Another advantage of the 
present approach, is that discretization oscillations are 
far smaller for local quantities appearing in E- lmv than 
for extensive quantities, such as (H) and ((H — (H)) 2 } 
appearing in the conventional approach to specific heats. 
Figure 3 shows the specific heat and entropy calculated 
with the above method, for the same parameters as in 
Fig. 1-2, with and without z-averaging. One sees that the 
discretization oscillations in the case of no z-averaging 
(n z = 1 curves) are drastically smaller than for the corre- 
sponding n z = 1 results from the conventional approach 
in Fig. 1. Including z-averaging makes the results of the 
new procedure indistinguishable from the Bethe ansatz 
calculations, as will be discussed in detail in Sec. IV-V. 



6 



2.0 
0.0 
"o-2.0 
m -4.0 
-6.0 
-8.0 



0.5- 



m 
it 

o o 



-0.5 



-(a) 



-| — i i i ii 



^imp 




.- E 

occ 




^docc 




hyb 






imp 



docc 
'hyb 



ill I I I I I 



ill I I I I I I I "I 

10 8 10" 6 10" 4 10" 2 10° 

k B T/D 



FIG. 4. (Color online) (a) The individual contributions 
E occ , Edocc, and -Ehyb to -Ej mp as a function of temperature 
(in units of Ao) for the symmetric model with parameters as in 
Fig. 1 (z-averaged with n z = 2). (b) Temperature derivatives 
of the above, yielding the relative contributions C occ , Cic.cc, 
and Chyb to the specific heat Ci mp . 
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FIG. 5. (Color online) (a) The individual contributions 
E occ , Edc.cc, and -Ehyb to -Ei mp as a function of temperature 
(in units of Ao) for the asymmetric model with parameters as 
in Fig. 1, but for an asymmetric level position Ed/Ao = — 1 
(z-averaged with n z = 2). (b) Temperature derivatives of 
the above, yielding the relative contributions C occ , Cdocc, and 
Chyb to the specific heat Ci mp . 



In Fig. 4(a), we show the different contributions -Eocc, 
^docc and -Ehyb /2 to the impurity internal energy for the 
symmetric Anderson model. Their temperature deriva- 
tives C occ ! Cdocc and Ch y b give the relative contribu- 
tions of these terms to the impurity specific heat C; mp 
and are shown in Fig. 4(b). Notice, that the Kondo in- 
duced peak in C mp at low temperatures results from a 
delicate balance of the hybridization (Ch y b) and Coulomb 
contributions (Cdocc), while the peak due to the resonant 
level at high temperatures is mainly due to the Coulomb 
term. The latter trend persists also for the asymmetric 
model, as shown in Fig. 5. Notice also that the gain in 
energy due to hybridization diminishes at high tempera- 
tures, reflecting the decoupling of the impurity from the 
conduction electrons in this limit. In general, however, 
the interaction of the impurity with the environment via 
the hybridization term provides an essential contribution 
at all non-zero hybridization strengths. 

(2) 

We now quantify the error in neg lecting dEr>(T)/dT 
in Eq. (16) for the calculation of impurity specific heats 
by, (a), comparing the result for Ci mp obtained within 
the new method with the Bcthc ansatz calculations, and, 

(2) 

(b) , explicitly calculating the contribution dE^ (T) jdT. 
Figure 6(a) shows the comparison to the Bethe ansatz 
calculation, where we also include the specific heat from 
the conventional approach. The relative deviation of 
the NRG calculations to the Bethe ansatz, shown in 
Fig. 6(b), is below 1% for all temperatures T < 0.01 = 
10 A . For T <C Tk, the relative error in C mp from the 
internal energy is 0.1% and 0.5% in the conventional ap- 
proach. The relative error exhibits remnants of the dis- 
cretization oscillations, which are not completely elimi- 
nated with z-averaging. Notice also that the errors in the 
two NRG calculations have the same error (relative to the 



Bethe ansatz) in the high temperature limit, T ^> Ao- 

(2) 

Hence, the latter error is not due to neglect of E^J in 
Eq. (15). Instead, it reflects, (a), the different high energy 
cut-off schemes in NRG and Bethe ansatz, and, (b), the 
finite size errors in the high energy excitation spectrum 
in NRG, since the latter stem from the shortest chains di- 
agonalized (typically m = 4 — 6) , which are also the ones 
most sensitive to the logarithmic discretization. The fact 
that the errors in both NRG calculations also correlate at 
lower temperatures (T < Ao) suggests that the neglect 

(2) 

of -Ej n( r in Eq. (15) is not the main source of error in cal- 
culating C mp (T). An explicit calculation that illustrates 

this is shown in Fig. 7. As stated above, the value of E-^J 
is of order A a /n, however, one clearly sees in Fig. 7(a) 

that E^2 has little temperature dependence (relative to 
the other contributions) for all temperatures extending 
up to the bandwidth D = 1. It's relative contribution to 
the impurity specific heat, shown in Fig. 7(b), for an en- 
ergy dependent A(oj), is negligible, typically contributing 
below 0.5%. 



IV. RESULTS FOR THE SYMMETRIC MODEL 

In this section we show results for the entropy and spe- 
cific heat of the Anderson model at the particle-hole sym- 
metric point Ed = —U/2. Results for zero magnetic field 
and increasing correlation strength U /Aq are presented 
in Sec. (IV A) and results for finite magnetic fields are 
given in Sec. (IV B). 

The symmetric Anderson model has been investigated 
in detail 1 and is well understood. For U/Aq ^> 1 and 
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FIG. 6. (Color online) (a) Comparison of specific heat, 
Ci m p(T), from the impurity internal energy (solid line) and 
conventional approach (dashed line) with the Bethe ansatz 
calculation (symbols) for the symmetric Anderson model with 
parameters as in Fig. I. NRG parameters also as in Fig. 1 
with n z = 2. (b) The relative deviation with respect to the 
Bethe ansatz result of the new (solid line) and conventional 
(dashed line) approaches. 



—Ed 3> Ao, a local spin S = 1/2 magnetic moment forms 
on the impurity. In this limit, the physics of the symmet- 
ric model at low temperatures T <§C min(|e,j + U\, \sd\,D) 
is that of the Kondo model 



Hk — H o + JS.s , 



(18) 



where, J is an antiferromagnetic exchange coupling be- 
tween the local spin S and the conduction electron spin- 
density s at the impurity site. The value of J is given by 
the Schrieffer- Wolff transformation J = 4V 2 /U. The 
low temperature properties (for U S> A ) are universal 
functions of T/Tk and B/Tk where we choose to define 
the Kondo scale from the Bethe ansatz result for the 
T = susceptibility x(0) via x(0) = (#Pb) 2 /4T k . For 
U > A , T K is given by 



T K = y/UA /2e 



-7rC//8Ao+7rA /2[/ 



(19) 



within corrections which are exponentially small in 
U/ttAq (see Ref. 1). For U = 0, the symmetric Anderson 
model reduces to a resonant level model and the relevant 
low temperature scale is then Aq. 



A. Zero magnetic field 

A comparison of the new approach with Bethe ansatz 
calculations is shown in Fig. 8 for the temperature de- 
pendence of the impurity specific heat and entropy for 
increasing values of the Coulomb interaction U / Ao . For 
U / A = 12, the Kondo induced peak in the specific heat 
at T p — aT^ with a ~ 0.29 is well separated from the 
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FIG. 7. (Color online) (a) The contribution E-^ to i?t mp 
as a function of temperature compared with E occ , Bdocc and 
Ehyb (in units of Ao) for the asymmetric model. Model pa- 
rameters: U — 12Aq, Ao = 0.001D, e^/Ao = —2 with a semi- 
elliptic hybridization function Im[A(tj)] = ~^d"\/ {D 2 — u> 2 ). 
A small A = 1.5 was used, which allows the spectral function 
entering E^ to be obtained without z-averaging. (b) The 
contribution Cdcriv = 8E^(T)/dT to Ci mp . The relative size 
of Cderiv to Ci mp lies between 0.2% and 0.5% for all temper- 
atures, except at temperatures approaching the bandwidth 
D = l. 



peak at T ps \sd\ due to the resonant level. With de- 
creasing U/Aq, the Kondo effect is suppressed and the 
Kondo induced peak in C{T) eventually merges with the 
peak due to the resonant level for U/Aq — > 0. Good 
agreement between the NRG and the exact Bethe ansatz 
calculations is seen for all values of U/Aq. 



B. Finite magnetic field 

At finite magnetic fields B > 0, the SU(2) spin sym- 
metry which we use in the NRG calculations, is broken. 
Therefore, in order to carry out calculations at finite 
magnetic field B > 0, preserving the numerical advan- 
tages of the full SU(2) symmetry, such as the increased 
number of states that can be retained, we obtained the 
finite field results by mapping the symmetric positive U 
Anderson model onto the negative-!/ Anderson model 
in the absence of a magnetic field but with local level 
given by £d = —U/2 — B/2 with U negative. 5 '' 58 This 
correspondence results from a particle-hole transforma- 
tion on the down spins only: — > d], d-j- — > df, and 



Ck^ with a particle-hole symmetric 



band e k = -e_ fe . 

Figure 9 shows the temperature dependence of 
Ci mp (T, B) for i3/Tx > 1 using our new approach and 
compared with Bethe ansatz calculations. The Kondo 
peak in the specific heat shifts to higher fields with in- 
creasing B and its position scales as B 2 /Tk for Tk <C 
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FIG. 8. (Color online) Temperature dependence of, (a), the 
impurity specific heat, Ci mp (T), and, (b), the impurity en- 
tropy, Si m p(T) / ln(2), for the symmetric Anderson model with 
Ao = 0.00 ID and increasing values of the Coulomb interac- 
tion: U I Ao = 4, 8, 12. Arrows in (a) indicated the Kondo 
scale 2k defined in Eq. 19. Symbols: new approach using 
NRG with A = 4 with an energy cut-off e c (A = 4) = 40, 
and ^-averaging [n z = 2, z = 1/4, 3/4]. Lines: corresponding 
Bethe ansatz calculations. 



FIG. 9. (Color online) Temperature dependence of the im- 
purity specific heat, Ci mp (T, B), for the symmetric Anderson 
model for U/Ao = 12, Ao = 0.00179 and increasing values 
of the magnetic field B/Tk > 1 where the Kondo scale T is 
defined in Eq. 19. Symbols: NRG calculations A = 4 with 
an energy cut-off e c (A = 4) = 40, and z-averaging [n z — 2, 
z — 1/4, 3/4]. Lines: Bethe ansatz calculations. Inset (a): 
T K -y(T, B) versus T/T K for several values of B/T K < 2, where 
7 (T,B) = C imp (T,B)/T. 



B <C £d- In contrast, the resonant level peak remains 
approximately fixed at T w Ed- As B approaches the 
value Ed, the two peaks merge into one peak at T sa Ed, 
with aproximately twice the height of the B — resonant 
level peak, and containing the whole entropy Si m p/&B = 
ln(4). The low field behaviour of Ci mp (T, B), also com- 
pared to Bethe ansatz calculations, is shown in Fig. 9(a) 
as Tk7(T, B) = C imp (T,B)/(T/T K ) versus T/T K for 
B/T K < 2. For T, B -> 0, 7 (T, B) -> 7 (0, 0) - 1/T K 
where 7(0, 0) is the linear coefficient of specific heat. This 
is strongly enhanced for U/A 3> 1 due to the exponen- 
tial decrease of Ik- A finite magnetic field of order 7k 
significantly suppresses the Kondo effect and results in 
smaller values of 7(0, B). As another check on the ac- 
curacy of our calculations, we estimate the Wilson ratio 
R w = 4tt 2 x(0)/37(0,0). This takes the value 2 in the 
Kondo regime of the symmetric Kondo model (i.e., for 
U ^> Ao) From the definition of Xk, we have that the 
susceptibility x(0) = 1/4T K , and from Fig. 9(a) we ex- 
tract 7(0,0) « 1.64/Tk, resulting in R w » 2.006, that 
is, a relative error in i?w below 1%. 
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FIG. 10. (Color online) Temperature dependence of Ci mp (T) 
for U/Ao ~ 12, Ao = 0.001D and local level positions £d/Ao 
ranging from Kondo (ed/Ao < —1), mixed valence (le^/Aol < 
1), and empty orbital (e,j/Ao > 1) regimes. Symbols: NRG 
calculations (new approach, ^-averaging and NRG parameters 
as in Fig. 1). Lines: Bethe ansatz calculations. 



V. RESULTS FOR THE ASYMMETRIC MODEL 

Figure 10 shows the impurity specific heat versus tem- 
perature for the asymmetric Anderson model, that is, 
for Ed > —U/2, calculated within the new approach. 
For comparison, we also show the corresponding Bethe 
ansatz calculations. One sees again excellent agreement 
at all temperatures between the two methods. Results for 



Ed < —U/2 are not shown, since these can be obtained 
from results for Ed > —U/2 by noting that the Ander- 
son model with parameters Ed, U, V transforms, under a 
particle-hole transformation applied to both spin species, 
to an Anderson model with parameters — (sd + U),U, V. 
This holds for a particle-hole symmetric constant density 
of states, the case considered here. 

The specific heat curves for the asymmetric model are 
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more complicated than those of the symmetric model. In 
the latter, the relevant excitations were the low temper- 
ature spin flip excitations, characterized by the Kondo 
scale Tk, and the excitations involving addition or re- 
moval of an electron from the resonant level, both char- 
acterized by an energy \ed\ = U/2. This accounts for the 
two peaks in the specific heat of the symmetric model: 
A high temperature peak at T ~ \sd\ and a low tempera- 
ture Kondo induced peak at T w Tk- For the asymmet- 
ric Anderson model, three types of excitation are pos- 
sible: Low-temperature spin flip excitations, associated 
with the Kondo scale T L = y/UAofa-*^ \^+u\/2UA 
of the asymmetric model, 1 and excitations associated 
with, (i), removing an electron from a singly occu- 
pied level (with energy scale \£d\) and (ii), removing 
an electron from a doubly occupied level (with energy 
scale | + U\). Thus, three peaks can be present in 
Ci m p(T): a Kondo induced peak at T w Tl, and two 
charge fluctuation induced peaks at T T\ = \e&\ 
and T fa T 2 = \ed + U\, respectively. In Fig. 10, the 
two high temperature peaks are seen in the mixed va- 
lence regime and partly also in the empty orbital regime 
(where the upper peak at T2 appears as a shoulder of 
the main peak at Ti). However, in the Kondo regime, 
the cases ed/A = —5, —3 with the choice U = 12A re- 
sult in Ti/A = 5,3 and T 2 /A = 7,9. In these cases, 
Ti and Ti are too close for separate peaks to be seen. 
In order to clarify this, we carried out calculations for 
U = 48A > A , and e d /A = -10,-8,-6,-4,-2 in 
the Kondo regime, for which Ti/A = 10,8,6,4,2 and 
T 2 /A = 38,40,42,44,46 > Ti/A are disparate scales. 
Figure 11 shows how the peaks at T ~ T\ and T « T 2 
evolve from the peak at T \ej\ — U/2 of the sym- 
metric model (dashed line in Fig. 11) on increasing e d 
above —U/2. Simultaneously, the Kondo peak in the spe- 
cific heat at Xl shifts to higher temperatures and eventu- 
ally merges with the peak at T\ when the mixed valence 
regime is reached (i.e., for Ed = — A ). Thereafter, only 
the high temperature peaks at T± and T2 are present. 
Notice also, that in the mixed valence regime T\ differs 
significantly from |e^|, a result of non-trivial renormal- 
izations present in the mixed valence regime, but absent 
in the empty orbital regime. 



VI. GENERALIZATION TO OTHER MODELS 



The approach of Sec. Ill can be straightforwardly gen- 
eralized to multiorbital and multichannel Anderson im- 
purity models with arbitrary local Coulomb interactions, 
as we briefly outline in Sec. VI A. In addition, in Sec. VI B 
we discuss it's application to dissipative two state sys- 
tems and the anisotropic Kondo model (AKM). 
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FIG. 11. (Color online) Temperature dependence of Ci mp (T) 
for (7/A = 48, A = 0.0001D and local level positions 
e d /A = -10, -6, -4, -2 (Kondo regime), e d /A = -1, 0, +1 
(mixed valence regime) and £d/Ao = +5 (empty orbital 
regime). NRG using the new approach (symbols) and con- 
ventional approach (solid lines) [A = 20, n z — 4,e c (A) = 
130]. Dashed line: resonant level peak in Ci mp at T/Aq « 
I Ed I / Ao = ?7/2Ao = 24 for the symmetric model (the Kondo 
induced peak at much lower T is not shown). The two high 
temperature peaks of the asymmetric model evolve from this 
peak when the asymmetry is finite. 



A. Multiorbital and multichannel Anderson models 

The multiorbital and multichannel Anderson impurity 
model is given by H = iJi mp + H + H int , where H inip = 
J2acr £ adi cr d ai7 + H C (U,U' , J), describes the impurity 
with a set of local levels having energies Sda, ot = l,...,g 
and Hq{JJ, U' , J) is the local Coulomb interaction involv- 
ing intra-orbital U, inter-orbital U' and a Hund's ex- 
change term J. The conduction electrons are described 
by H = J^kacr e k a c\ a<J c kaiJ where e ka is the kinetic 
energy of electrons in band a. These bands hybridize 
with hybridization strengths V a ,a = 1, ... ,g to the lo- 
cal levels via H int = J^kaa V a {cl acr d aa + S aa c kac! ). Let 
A a {ui) — J2k Va /( UJ ~ e ka) denote the hybridization func- 
tions characterizing H lnt . Proceeding as in Sec. Ill, we 
write the impurity internal energy as -Eimp = Btotal ~ Eq 
where -Etotai = (H) is the total energy and E = (H ) = 
J2kaa e ka(cl acr c kacr )o is the energy of the non-interacting 
conduction electrons in the absence of the impurity. The 
latter is given by E = ]T Q(T / def(e)eN a (e), where /(e) is 
the Fermi function and N a (e) — 5{e—e ka ) is the non- 
interacting conduction electron density of states per spin 
for band a. Stotai is a sum of local occupation number 
contributions -B cc = ^2 ar7 £ a( n atr) and local Coulomb 
terms Ec — (Hc(U, U', J)} and two further terms involv- 
ing the interacting band E cond = J^kaa e ka{c\ aa c ka(J ) 
and the hybridization energy £; hyb = J2aa V a (<ft aer foao + 
H.c.) where V a f 0aa = J2k c ka<r'- 

-Etotai = E occ + Eq + E COD( i + E'hyb- (20) 
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We evaluate the latter two contributions as in Sec. Ill, 
finding 



E hyh = f duf(u)Im[G da(T (Lj)A a (<j)} , (21) 

and -Econd = A + E- mt , where 



Ant = lY,J <^f(") h 



d 

G daa (uj) — (wA a (w)) 



, p(2) 
Ant + Ant > 



Ant =~J2 f dLof(oj)Im{G daa (Lj)A a (uj)} 

rurr 



<9w 



(23) 



(24) 



and Gdaa{u) is the retarded Green function for local level 
a. Combining E^ with £"hyb gives for the impurity in- 
ternal energy 



A 



imp — E occ + Eq + ^ Aiyb + Ant ) 



(25) 



where, as before, all contributions except the last one are 
evaluated as local static correlation functions. For rea- 
sons discussed in Sec. Ill, the temperature dependence of 
the last term is negligible in many cases and the impurity 
specific heat can be calculated to high accuracy via 

1 dE hyh 



dE n 



dEr 



imp " + fW" 2~dT 



dE 



ionic 

dT 



1 *9Aiyb 

2 dT : 



(26) 
(27) 



where E-„ 



B. 



(H 



imp/ 



Dissipative two state systems and the 
anisotropic Kondo model 



The method of Sec. Ill can be applied to bosonic mod- 
els such as the dissipative two state system, 4 ' and for 
Ohmic dissipation, one can further relate the results to 
the AKM and related models (for example, a two-level 
system in a metallic environment " ). Dissipative two 
state systems are of interest in many contexts, including 
the description of qubits coupled to their environment. 

The Hamiltonian of the dissipative two state system 
is given by H = H§ + Hb + Hi. The first term Hg = 
— \Aq<j x + \t&z describes a two-level system with bias 
splitting e and tunneling amplitude A , and a i — xvz are 

Pauli spin matrices. Hb = X)j ^>i( a t a i + 1/2) is the envi- 
ronment and consists of an infinite set of harmonic oscil- 
lators (i = 1,2,..., oo) with a.i(aj) the annihilation (cre- 
ation) operators for a harmonic oscillator of frequency 
u>i and < u>i < w c , where uj c is an upper cut-off fre- 
quency. The non-interacting density of states of the 
environment is denoted by g(u>i) = J+ <K W — w «) an d 



is finite in the interval [0,w c ] and zero otherwise. Fi- 
nally, H\ = ha z Aj(aj + aj) describes the coupling 
of the two-state system co-ordinate a z to the oscilla- 
tors, with Xi denoting the coupling strength to oscilla- 
tor i. The function T(uj + id) = 2iOi/2) 2 /(w - w l + 
iS) = J dw'(A(u/)/2) 2 g(uj')/(uj — u' + iS) characterizes 
the system-environment interaction. The Ohmic two 
state system, specified by a spectral function J{uj) = 
— — ImT(w + iS) ~ auj for cj — > 0, where a is the dimcn- 
sionless dissipation strength, is equivalent to the AKM 

-S-s+)+J»S z s% + BS z , 



(22) H = T,ka ^c\ a c ka 



where Jj_ (Jn is the transverse (longitudinal) part of the 
Kondo exchange interaction and B is a local magnetic 
field. The correspondence is given by pj± = — Aq /uj c 
and a = (1 + 26 /ir) 2 where 5 = arctan(— 7rpJii/4) and 
p is the density of states of the conduction electrons in 
the AKM. 4,5,60-62 The low energy scale of the Ohmic 
two state system is the renormalized tunneling ampli- 
tude A r given by A r /u; c = (Ao/w c ) 1 ^ 1 ~ Q '' and corre- 
sponds to the low energy Kondo scale Xk of the AKM. 
Special care is needed to obtain results for the Ohmic 
two state system from the AKM in the vicinity of the 
singular point a — > l~, since this corresponds to J|| — > 
but with the condition < J± < Jn , that is, in terms of 
parameters of the Ohmic two state system one requires 
Ao/w c <C 1 — a <C 1 in order to investigate the vicinity 
of a = 1 within the AKM. 5 



The specific heat, G\ 



imp 



zp — 

1) is 



8E imp /dT, of the Ohmic 
two-state system is defined via an impurity internal en- 
ergy E imp = Aotal - #0, where Aotai = (H) = (H s ) + 
(H B ) + (Hi) and E = (H B )o = E^M^O + A 
Jq C doj uj n(uj) g(uj) + A P where n(uj) — l/(e^ u 
the Bose distribution function and the zero point en- 
ergy E zp can be dropped, as it cancels in the difference 
(Hb) — Eq = Eb — Eq appearing in E lmp . Evaluating 
Eb—Eq and Ei = (H) following the approach in Sec. Ill, 
we find 

d 



Eb — Eq = 



E 



(i) _ 



A 



(2) 



dujn(uj)lm 



X: 



+ id)-—- (ojT(u 
ouj 



(i) 



E 



(2) 



iS)) 
(28) 

dwn(oj)Im[x zz (u) + id)T(w + i5)] (29) 

dT(uj + iS) 



dum{uj)lTii 



Xzz{u + iS)uj- 



(30) 



and 



E, 



dwn(w)Im [x zz (u; + i5)T(uj + iS)] , (31) 

(32) 

where Xzz{w + iS) = ((<J z ]Vz))u+i8 is the longitudinal 
retarded dynamic susceptibility and T(co + iS), charac- 
terizing the system-environment interaction, was defined 
above. Noting that exactly cancels E\ in the impu- 
rity internal energy, we find 

1 



E u 



-A ((T X ) 



E 



(2) 



(33) 
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that is, E imp — Eg + E^' . The term E^ gives a non- 
negligible contribution to the impurity internal energy. 
For example, in the Ohmic case with spectral function 
J(oj) = — ^lmT(ui + i5) ~ olui we have u>dJ(uj)/duj ~ aw 

(2) 

at low frequencies, so E^ provides a contribution pro- 
portional to a. By carrying out specific heat calculations 
on the AKM, we find numerically that the impurity spe- 
cific heat is consistent with setting E^ = haAo(cr x ) +A, 
with A being a weakly temperature dependent term, and 
negligible for calculating the specific heat, except in the 
limit a — > 1~ . The latter limit is difficult to treat nu- 
merically because of the vanishing low energy scale A r 
for a — > l - (e.g., for Aq/u c = 0.01 and a = 0.9 we have 
A r /uj c = 10 -20 ). Hence, except in this extreme limit, and 
as we show below by comparing with exact results, the 
impurity specific heat can be obtained accurately from 
Cimp = i5 §p by using 

E imp &-^A Q (l-a)(a x ) + ^e(a z ). (34) 

Figure 12 shows results obtained in this way for 
Cimp(r) / (kftT I A r ) compared to Bethe ansatz calcula- 
tions for the AKM ,,,! for a range of dissipation strengths. 
These results recover the known results for asymptoti- 
cally high and low temperatures. 64 In common with spe- 
cific heats of other correlated electron systems as a func- 
tion of interaction strength," we observe a crossing point 
in C{T)/T (here, at k B T/A r 0.67). On decreasing 
the dissipation strength from strong (a > 1/2) to weak 
values (a < 1/2) the T 3 coefficient of the specific heat 
changes sign for a < 1/3 resulting in the appearance of 
a finite temperature peak in C(T)/T. This is shown in 
more detail in Fig. 13. It signifies the development of a 
gap ~ A in the spectrum as a — > 0. For a = one 
eventually recovers the Schottky specific heat for a non- 
interacting two level system. The expression (34) for the 
Ohmic two system is also the impurity internal energy 
of the equivalent AKM (indeed, the NRG results that 
we showed were for this model). The correspondence of 
model parameters was given above and the operators o~ x 
and a z are identified, under bosonization, 4 ' 5,60 ' 62 with the 
spin-flip operator S + Sq + S~s~q and the local S z in the 
AKM, respectively. The zero temperature expectation 
values (a x ) and (a z ) (and the associated entanglement 
entropy of the qubit) have been studied previously as a 
function of dissipation strength and finite bias. ' 

(2) 

We expect that the term Ex ' is non-negligible also for 
generic spectral functions J(lS) ~ w s and certainly for the 
sub-Ohmic case s < 1. Recent results for the local spin 
dynamics of the sub-Ohmic spin boson model could 
shed light on this. 

The result (34) shows that a significant contribution to 
the impurity internal energy and specific heat arises from 
the (interacting) bath contribution E B , which remains 
finite for arbitrarily small a. Thus, while a definition 
of the internal energy of the system via E$ = (Hg) and 
the specific heat via Cs = dE§/dT, might seem reason- 
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FIG. 12. (Color online) Specific heat, G mp (T)/fc B T/A rj 
of the Ohmic two state system as a function of reduced 
temperature ksT/A T , for a range of dissipation strengths 
= 1/5,1/4,1/3,1/2,2/3,3/4,4/5. Symbols: NRG results 
in new approach. Lines: Bethe ansatz results. The renormal- 
ized tunneling amplitude A r from the Bethe ansatz is used. 
The vertical arrow indicates the approximate crossing point 
at k B T/A T w 0.67. Model parameters: A /oj c = 0.005. NRG 
parameters: A = 10, n z — 4 retaining 860 states per NRG 
iteration. 

able for a small quantum system weakly coupled to an 
infinite bath, such a definition yields, in general, a spe- 
cific heat Cs which differs from Ci mp . <,9 ~' 2 One system 
for which the two definitions agree is the harmonic oscil- 
lator coupled Ohmically to an infinite bath of harmonic 
oscillators. (>l> This result, however, represents a special 
case, and, moreover, is sensitive to details of the cut-off 
scheme used for the spectral function J{oS) (see Ref. 09 
and 72). The use of E- lmv and Ci mp as definitions for the 
system internal energy and specific heat in the context 
of open quantum systems ' also provides an unambigu- 
ous prescription for their measurement in terms of two 
separate measurements, ' one for H and one for Hq. 

We note also that the impurity specific heat Ci mp (T) = 
C(T) — Cq(T) need not be positive at all temperatures 
and only the positivity of C(T) and Cq(T) in Eqs. (1) 
and (2) is guaranteed by thermodynamic stability of the 
equilibrium systems described by H and H (see Ref. 75). 
Examples of systems where the difference, Ci mp (T), may 
be negative in some temperature range, include quan- 
tum impurities exhibiting a flow between a stable and 
an unstable fixed point, and magnetic impurities in 
superconductors . 

VII. DISCUSSION AND CONCLUSIONS 

In this paper, we introduced a new approach to the cal- 
culation of impurity internal energies and specific heats 
of quantum impurity models within the NRG method. 
For general Anderson impurity models, the impurity con- 
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FIG. 13. (Color online) Specific heat, C imp (T)/(ak B T/A T ), 
of the Ohmic two state system as a function of reduced tem- 
perature fcsT/A r , for a range of dissipation strengths a < 
1/2. Symbols: NRG results in new approach. Lines: Bethe 
ansatz results. In the low temperature Fermi liquid regime, 
T < A r , we have C imp (T)/(afe B T/A r ) = 7 + /3(T/A r ) 2 , 
with 7 = 7r/3 and the T 3 coefficient in C(T) changes sign 
for a < 1/3 (see Ref. 62). Model parameters: Aq/uj c = 0.005. 
NRG parameters: A = 10, n z — 4 retaining 860 states per 
NRG iteration. 



tribution to the internal energy was expressed in terms 
of local quantities and the main contribution to the im- 
purity specific heat was shown to arise from local static 
correlation functions. For this class of models, the im- 
purity specific heat can be obtained essentially exactly 



as Cimp(^) — gr" 1 2 ' wnere -^ionic — (-Simp) 

and .Ehyb is the hybridization energy. A comparison with 
exact Bethe ansatz calculations showed that the results 
for specific heats of the Anderson impurity model are re- 
covered accurately over the whole temperature and mag- 
netic field range. The new method has several advantages 
over the conventional approach to specific heats within 
the NRG, namely, (i), only diagonalization of the total 
system is required, (ii), only local quantities are required, 
and, (hi), discretization oscillations at large A are signif- 
icantly smaller than in the conventional approach. 

For the dissipative two state system we obtain the spe- 
cific heat as C; mp (T) = dE '^ p — + 9 gj, , where 
Eg = (Hs) is analogous to iconic in the Anderson model, 

(2) 

and -Eg is a contribution to the energy of the system 
arising from the interaction with the bath. It depends 
on the local dynamical susceptibility and the type of cou- 
pling to the environment. For the Ohmic case, we used 
the equivalence of the Ohmic two state system to the 

(2) l 

AKM to show numerically that E^' = ^aA (a x ) + A 
with A having a negligible temperature dependence, ex- 
cept in the extreme limit a — > l~. Comparison with 
exact Bethe ansatz calculations on the AKM confirmed 
the above. 

The approach described in this paper applies to en- 



ergy dependent hybridizations also, see Fig. 7, so, in- 

(2) 

elusion of the term E> n £ in Eq. (15), could prove useful 
in applications to quantum impurities with a pseudogap 
density of states.' 5 ''' 8 It may also be applied within other 
methods for solving quantum impurity models, for exam- 
ple, within continuous time or Hirsch-Fey quantum 
Monte Carlo techniques or exact diagonalization methods 
(for a recent review see Ref. 81 and references therein). 
Local static correlation functions, such as the double oc- 
cupancy, required for E- lmp , arc readily extracted within 
these approaches. 82 

Within a DMFT treatment of correlated lattice 
models, " the hybridization function A acquires an im- 
portant temperature and frequency dependence A(w) —¥ 

(2) 

A(u>,T). The latter enters explicitly in the term E> nt ', 
whose inclusion could offer an approach to the calula- 
tion of specific heats of correlated lattice models. The 
thermodynamic potential of the latter 8 ' is a sum of 
two parts, one depending on the local self-energy, which 
is the central quantity calculated in DMFT, and an- 
other equal to the thermodynamic potential, f2i mp = 



Es 



TS; r 



imp 

The 



imp ~ -l ^imp j of the effective impurity model 
latter can be obtained from Ei mp (T), via C lmp (T) and 

S imp (T) = fo dT' Cim ^ T,) . The impurity internal en- 
ergy, expressed in terms of local dynamical quantities as 
in Ref. 45, has recently been used in a DMFT solution of 
the Hubbard model within a variational generalization 
of the local moment approach. 8 '' 

In the future, it may be interesting, especially in the 
context of qubits or nanodevices, to consider the time 
dependence of the impurity internal energy subject to an 
initial state preparation, for example, within techniques 
such as time-dependent density matrix renormalization 
group - or time-dependent NRG.' 5 ' 3 ' 89 ' 90 
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Appendix A: Band contribution to impurity internal 
energy 

The expression (11) for the conduction band contribu- 
tion to the impurity internal energy requires evaluation 
of the integral 



I(uj) = j de 



eV 2 N(e) 
(u> - e + i5f 



(Al) 



We assume a density of states N(ui) vanishing at the band 
edges at u> = ±D. The hybridization function A(w) = 



13 



J2k V 2 /(^ - £fe + iS) = A fl (w) + iAj(w) where Ai(w) 
— 7rA^(w)F 2 . With these definitions, we have 



+D 



deeA/(e) 



1 eAj(e) , +D 

i r 1 ' . i 



7r w — e + « 

f-D 

de 

D 



7T 



<9e (w 



i<5 <9e 



(eA^e)). (A2) 



The first term vanishes since Aj(±D) = for regu- 
lar (e.g., 3D) densities of states (and will otherwise re- 
sult in contributions with negligible temperature depen- 
dence). The second term can be evaluated by noting 
that A(w + iS) satisfies the causal properties of retarded 
Green functions and by using the following properties of 
principle value (P.V.) integrals: if P.V.[/(a;)] = g(y) then 
P.V.[/'(z)] = g'(y) andP.V.[x/(x)] = yg(y)+± J dxf(x). 
The final result is 



I(w) 



( W A(w)) 



(A3) 



Appendix B: Numerical solution of the 
Thermodynamic Bethe Ansatz equations 



In this Appendix, we summarize the thermodynamic 
Bethe ansatz (TBA) equations for the Anderson model, 
which were derived by Okiji and Kawakami 42 ' 91 ' 92 and 
Tsvelick, Filyov, and Wiegmann, 40 ' 41 ' 9! ' 94 and provide 
details of their numerical solution. 46-49 ' 62 ' 95 The numer- 
ical procedure described applies to both the symmetric 
and asymmetric Anderson models and in the presence of 
a finite magnetic field and was used to obtain the results 
presented in this paper. 



1. Thermodynamic Bethe Ansatz Equations 



The thermodynamic Bethe ansatz (TBA) produces an 
infinite set of coupled integral equations for the func- 
tions e(fc), n' n (A) and K n (A), n — 1,2, ... , describing the 
charge and spin excitations of the system (Tsvelick and 



Wiegmann 40 ): 

/oo 
s(g(k)-A)-]n(f(K 1 (A))dA = 
-OO 

eo(fc) — T ■ f°° s(g{k) - A) ■ ln(/«(A))dA 

J — OO 

(Bla) 

Kn (A) +T-(s* (ln(/(/c„ +1 )) + ln(/( Kn _ 1 ))))(A) = 

Sn,i T - [°° s(g(k) - A) • ln(/(-e(fe))) • g'(k) dk 

(Bib) 

<(A) + T • (s * (ln(/« +1 )) + ln(/«_ x ))))(A) = 

/OO 
s(g(k) - A) -ln(f(e(k))) ■ g'(k)dk 
-OO 

(Blc) 

where 



9{k) = 
f(k) = 



{k-e d ~\U) 2 

2ru 



l + e fe / T ' R{ ~ X) n 



s(A) = 
1 



1 



2 cosh(7rA) 
cos(o;a;) 



1 +e u 



du 



e (k) = k-e d --U + 



R{g{k) - g(p)) -p-g'{p)dp 



g'{k) denotes the first derivative of g(k) with respect to 
k. * is the convolution of two functions, kq and k' equal 
— oo. For 7i — y oo the functions approach the constant 
values, 

lim K n = n ■ H, lim k' = n ■ (2ed + U), (B2) 

n— foo n— >oo 

where H is a uniform magnetic field and 2ed + U mea- 
sures the deviation from the symmetric point at Sd — 
— U/2. The impurity contribution to the specific heat, 
Cimpj may be calculated from the the impurity contribu- 
tion to the thermodynamic potential, r2; mp , via Ci mp = 
-Td 2 fl imp /dT 2 : where 

/OO 
p Q (k) • ln(/(-e(fc)))dfc 
OO 

/OO 
fro(A)-m(/ti(A))dA + -Eo (B3) 
-OO 

The functions po & n d <7o are given by: 

/oo 
s(A - g(k)) ■ A(k)dk 
-OO 

/OO 
R(g(k) - g(p)) ■ A(p)d P , 
-oo 

where A(fc) = ^^ V 2 + ^ k _ ed ^ ■ Eq is the ground state 

energy of the symmetric Anderson model. Note two 
changes with respect to the earlier Ref . I : a sign change 
in equation Blc (as in Wiegmann and Tsvelick ) and a 
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FIG. 14. (Color online) The figure shows a set of £ n for 
the symmetric case (e<2 + U/2 — 0) zoomed to range of A = 
— 20 . . . 20. The functions become smoother with higher n due 
to the convolution with s(x). 



factor 2 in the boundary value for n' n in equation B2 (as 
in Okiji and Kawakami 42 ). 

For the calculations we use a transformation of n n and 
n' n to new functions £ n and £f n similar to that used in pre- 
vious works. 46,4 '' 62 After substituting = ln(l + e K ™/ T ) 
and & = ln(l + e</ T ) we obtain the following coupled 
equations: 



6(A) 
f»(A) 
£(A) 
C(A) 

h(A) 
I'M) 



ln(l + exp((s*(6+/i))(A))) (B4a) 

ln(l + cxp(( S * (e„-i + e»+i))(A))) (B4b) 

ln(l + exp(( S *(^+70)(A))) (B4c) 

ln(l + cxp(( S * + C+i))(A))) (B4d) 

S ( fl (fc)-A)-ln(/(-e(fc)))- 5 '(fc)dfe 

OO 

(B4e) 

X) 

S ( 3 (fc)-A)-ln(/(e(fc))). 5 '(fc)dfc 

OO 

(B4f) 

X> 

»( fl (fc)-A)-(6(A)-ri(A))dA (B4g) 



!(*) = 

e(k) = e (k) + T-I(k) 



(B4h) 



for £jv and £' N the corresponding truncation functions are 
calculated by: 



In 



((cOSh(y)- ^ 



^l + siim 2 (|).[2 + e«"-i]) 2 ) 
In ^( cosh( 



(B5a) 



,2e d + U , 
2 ' 



2 + c 5 «'-i- 



l + sinh 2 (?£±±^).[2 + e«« 



- l ]) 2 ) (B5b) 



As a further check, and to ensure the correct be- 
haviour at the boundaries, the TBA integral equations 
were explicitly solved in the limits of A, k — > ±oo. As 
the functions are smooth in this limit one can assume 
that s(x) -> S(x)/2 and lim e (fc) = 2(k - e d - U/2), 

k— foo 

lim e (k) = 0. This leads to the following set of cou- 

k — ^ — oo 

pled algebraic equations: 
lim 



A— >■— oo 




l-exp(ifc))) 






= ln(l - 


(B6a) 




= ln(l - 


1- exp(-(^ Tl _i 


+e„+i))) (B6b) 




= ln(l - 




(B6c) 




= ln(l - 


^exp(-(^_! 


+ (B6d) 


lim 








A— >-oo 




^ ex p(^(6 - 


ln(l+cxp(^i))))) 




= ln(l - 






1- exp(-(^ T1 _! 


(B6c) 




= ln(l - 


+ £ n+1 ))) (B6f) 




= 




(B6g) 




= ln(l - 


+ C+i))) (B6h) 



The truncation constants £jv and are calculated as 
in equation B5. The boundary values where calculated 
by iteration using a modification of the Powell hybrid 
method. 



2. Truncation 

For calculational purposes the equations and 
have to be truncated at some finite value n = N and 
n' = N' . One has to calculate the functions at the trun- 
cation with care, to avoid wrong results at the boundaries 
A — » ±oo. We use the truncation scheme of Takahashi 
and Shiroishi. It is assumed that the function s(x) can 
be approximated by 5{x)/2 for large n or nf. This is 
justified as the functions become smoother in this region 
(see figure 14). Rewritten for the Anderson Model and 



3. Numerical Details 

For the calculations, a logarithmic grid was used that 
is centred around Ed + U/2. The TBA equations were 
solved by iteration. The initial values of £ n and / were 
chosen to fit a tanh-function with boundary values given 
by the correct boundary values of and £' n , obtained 
as described above. The integrations were carried out 
using adaptive routines with the integrands being repre- 
sented by splines of smooth functions only (see below). 
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FIG. 15. (Color online) Comparison between I[ and for 
500 iterations and T = 10 _4 Tk- Parameters were chosen 
to be the same as in figure 1. For very low temperatures 
£i (circles, left y-axis) exhibits an exponential drop beyond a 
certain rapidity Ao (~ —4 for the case shown) which is difficult 
to capture with a fixed grid. This problem can be overcome 
by using the smooth function I[ (squares, right y-axis) to 
calculate via Eq. (B4c). 



A smoother convergence of the iteration procedure is ob- 
tained by using 10% of the old iteration values in each 
step. To represent only smooth functions as splines, £i 
and are not interpolated, but instead the s * £2 and 
s * £ 2 respectively. The values of £1 and (;[ are then cal- 
culated from these convolutions and from 1\ and /{ using 
Eq. (B4a) and Eq. (B4c) . This avoids numerical problems 
due to the exponential drop to zero of Ci beyond a certain 
rapidity Ao. See Figure 15 for a comparison between the 
behavior of !;[ and I[. N = N' = 20 functions were used 
and iterated 500 times for the figures in this section (and 
2000 times for the results in the paper) . The growth-rate 
of the grid was 1.05 and it consisted of 801 points. The 
mid 400 values lie in a range of [—40,40]. After a cer- 
tain temperature dependent cut-off (±40 ± 40 • T/T ) the 
boundary values were used instead of being calculated to 
ensure numerical stability. The thermodynamic poten- 
tial was calculated in a range of To • 10~ 3 to To ■ 10 6 on a 
logarithmic mesh (factor 2 1 / 8 as step width) where To is 
defined as T = ^C/r/2-exp(-|^ + Kondo temper- 
ature for the symmetric case. It is related to the magnetic 

susceptibility at zero temperature Xi mp (T = 0) = ^ B j o 
(see Hewson in Ref. I p. 165, and Kawakami and Okiji 
in Ref. 96). 
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